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Abstract 



Gap solitons near a band edge of a spatially periodic nonlinear PDE can be formally approximated 
by solutions of Coupled Mode Equations (CMEs). Here we study this approximation for the case of 
the 2D Periodic Nonlinear Schrodinger / Gross-Pitaevskii Equation with a non-separable potential 
of finite contrast. We show that unlike in the case of separable potentials [T. Dohnal, D. Pelinovsky, 
and G. Schneider, J. Nonlin. Sci. 19, 95-131 (2009)] the CME derivation has to be carried out 
in Bloch rather than physical coordinates. Using the Lyapunov-Schmidt reduction we then give a 
rigorous justification of the CMEs as an asymptotic model for reversible non-degenerate gap solitons 
and provide if" estimates for this approximation. The results are confirmed by numerical examples 
including some new families of CMEs and gap solitons absent for separable potentials. 

1 Introduction 

Coherent structures, like gap solitons, in nonlinear periodic wave propagation problems are important 
both theoretically and in applications. Typical examples include optical waves in photonic crystals 
and matter waves in Bose-Einstein condensates loaded onto optical lattices. A standard model in 
these contexts is the Nonlinear Schrodinger/Gross-Pitaevskii equation with a periodic potential, which 
applies in Kerr-nonlinear photonic crystals |37|, \TE[ [2^1 [TB] as well as in Bose-Einstein condensates 
loaded onto an optical lattice [201 EH • Here we consider the case of two spatial dimensions and 
without loss of generality take the potential 27r-periodic in both directions and, hence, consider 



iEt = -AE + V{x)E + cr\E\^E, V{xi+2tt, X2) = F(xi, X2+27r) = V{x), a; G t G M, (1.1) 

with E = E{x, t) eC,a = ±l and V G i^j™ (M^), > 1, m G M. 

We are interested in stationary gap solitons (GSs) E{x,t) = (/>(a;)e~'^*. Thus 4> solves 
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where soliton is understood in the sense of a sohtary wave, which means that |</'(x)| — )• exponentiahy 
as |x| — )• 0. Necessarily, then uj has to he in a gap of the essential spectrum of the operator L := 
—A + V{x), hence the name "gap soliton." From the phenomenological and experimental point of view 
multidimensional GSs have been widely studied in the context of both photonic crystals [U [25l [151 [13 
[T2| [T6] and Bose-Einstein condensates [H [261 IZl • 

The essential spectrum of L is given by the so called band structure, and for our analysis we choose 
cj close to a band edge, i.e., u = + e^fi, < e <C 1, where is an edge of a band gap and O 
has a sign chosen so that ui lies inside the gap. Using a multiple scales expansion one may formally 
derive coupled mode equations (CMEs) to approximate envelopes of the gap solitons near gap edges. 
CMEs are a constant coefficient problem formulated in slowly varying variables. They are, therefore, 
typically more amiable to analysis and also cheaper for numerical approximations compared to the 



original system (1.2). The multiple scales approach has been used both for the Gross-Pitaevskii and 
Maxwell equations with infinitesimal, i.e. 0{e), contrast in the periodicity V{x) [21 [6l [Ml [3l [5l [12] 
as well as with finite contrast |111 l33l [T3] . The main difference in the asymptotic approximation of 
the two cases is that for infinitesimal contrasts the expansion modes are Fourier waves while for finite 
contrast they are Bloch waves. However, in dimension two and higher sufficiently large (finite) contrast 
is necessary to generate band gaps due to overlapping of bands in the corresponding homogeneous 
medium. The only exception is the semi-infinite gap of the band structure of L of the Gross-Pitaevskii 
equation. As a result, gap solitons in finite gaps of the Gross-Pitaevskii equation and in any gap of 
Maxwell systems in dimensions two and higher can only exist for finite contrast structures. 

Localized solutions of CMEs formally yield gap solitons of the original system. However, the formal 
derivation of the CMEs, discarding some error at higher order in e, does not imply that all localized 
solutions of the CMEs yield gap solitons. For this we need to estimate the error in some function 
space and to show the persistence of the CME solitons under perturbation of the CMEs. A famous 
result concerning non-persistence is the non-existence of breathers in perturbations of the sine-Gordon 
equation, e.g., [10]. On the other hand, GSs are known to exist in every band gap of L, see, e.g. [361127]. 
The proofs, however, are based on variational methods and do not relate GSs to solutions of the CMEs. 



A rigorous justification of the CMEs has been given for (1.2) in ID in [28J, and in 2D in [13| . but 
only for the case of a separable potential 



V{xi,X2) = Wi{xi) + W2{x2). 



Here we transfer this result to not necessarily separable potentials, where we need some minimal 
smoothness, namely V £ H^^(M?),m > 1. As an example we choose 



V{XI,X2) = l + {7]-l)W{xi)W{x2), W{S) 



tanh I 7 



) + tanh ( 7 



Svr 



, (1.3) 



which represents a square geometry with smoothed-out edges (Fig. [T]). We choose the contrast r]mV 
so that two finite gaps appear in the band structure of the corresponding linear eigenvalue problem. 
One main difference between the separable and non-separable case lies in the fact that for the non- 
separable case band edges may be attained at wavenumbers not within the set of vertices of the first 
irreducible Brillouin zone. Then the CME derivation and justification is impossible to carry out in 
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V(x^,X2) 




Figure 1: The periodic potential V in (1.3) over the Wigner-Seitz cell. 



physical variables and has to be performed in Bloch variables. This case occurs at least at one band 



edge of the potential (1.3), and the presented CMEs corresponding to this edge have, to our knowledge, 
not been studied before. Similarly, the GSs which we show to bifurcate from this edge are new. 

In f|2] we discuss in detail the band structure for (1.3) and the associated Bloch eigenfunctions. 



together with their symmetries. Then in ^we first give the formal derivation of the CME in physical 
space, reporting a failure in one case where the band edge is attained simultaneously at four wave 
numbers outside the set of vertices of the first Brillouin zone, and present a general CME derivation in 
Bloch variables. The existence of gap solitons is proved in phased on the existence of special (namely 
reversible and non-degenerate, see below) localized solutions of the CMEs, in the following sense. 
Function spaces and notation. For m G N, the Sobolev spaces H"^(M?) are classically defined as 

if™(]R2) := {u G L^{R'i) : d^u £ L^iR"^) for \a\ < m}, with norm \\u\\hs = (E|a|<m ll^x^^llia) ^ ^ 
where d^u denotes the distributional derivative, see, e.g., [Ij. Then, for s = m G N, Fourier transform 

^{k) := {TcP){k) := f </.(x)e-'^-^-dx, 4>{x) = {T^^ m) ■-= I me'"-^ dk, (1.4) 

is an isomorphism from H^(M?) to 

l2(m2) := e L\R') : Uh^^ := ||(1 + |fc|)>||i2 < oo}, i.e. CMl2 < UU^ < C2m^. (1.5) 

From the applied point of view we could restrict to integer s. However, since our analysis is strongly 
based on Fourier transform, it is conceptually cleaner to use a definition of Sobolev spaces based on 
with arbitrary s > 0. Thus, henceforth we use H^(M?) := J^^^L^(M^) as definition for < s G M. This 
also gives a very simple proof of the Sobolev embedding theorem < C||(/)||/fs for k < s — 1, see 

Lemma 14.21 below. 

Main result. Let s > 2 and V G H^^^ ^^^{M.'^), 6 > 0, where [s] is the smallest integer larger than or 
equal to s. Let A = (^i, . . . , A^) be a reversible non-degenerate localized solution of the CMEs with 
A G [H''(M?)]^ for ah g > 0. Then for co = + e'^Q with sufficiently small there exists a GS (pes 
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for (1.2), such that (pes G H''(R'^), and (/>g5 can be approximated by 

N 

= e^Aj{ex)un,{k^''^;x), (1.6) 
i=i 

where Unj{k^-^^; x) G i7j[^*J''^^^'^(M^) are the pertinent Bloch waves, j = 1, . . . , N. In detail, we prove 

UGS-e<P^'^^\HHR^)<Ce^/^ (1.7) 

where the estimate can be improved in special cases, see below. 

Note that ||e(^(°) 1 1^00(^2) = 0{e) but || 11^2(^2) = 0(1) such that the error in il.7n is indeed 
smaller than the approximation. The proof is based on a Lyapunov-Schmidt reduction and analysis 
of suitable extended CMEs. In ^we give some numerical illustrations and verify convergence of the 
asymptotic coupled mode approximation. 

Remark 1.1 The (apparent) lack of an estimate 

Ugs - e0(°)||Loo(K2) < Cei+^ with /? > IS a 
disadvantage of our analysis. It is due to the fact that we work in L^-spaces in Fourier resp. Bloch 
variables, while a direct L°° estimate in physical variables would require working in L^-spaces in Fourier 
resp. Bloch variables. This is not possible due to a technical obstacle, see \13\ §8]. On the other hand, 
Hilbert spaces are also more natural spaces to work in since they allow direct transition from physical 
to Bloch variables and back. Moreover, localization in x in the sense of decay to for |x| — )• oo follows 
directly in spaces of integrable functions. Note also that based on the formal asymptotics, instead of 
£4/5 Qj^g pg^j^ expect the convergence rate in H'^ in (1.7) which is the approximate rate observed in 



our numerical examples. Finally, in Remark |4.13| we explain how the long wave modulational form of 
the formal asymptotics allows to obtain an 0{e^^^) convergence of the error in from the 0{e^) 
convergence in H^. However, a completely rigorous calculation is lengthy and therefore here we content 



ourselves with (1.7). J 



Remark 1.2 Time-dependent CMEs have been justified in ID for infinitesimal |19| I32j and finite [8] 
contrast, and in 2D for finite contrast under the condition of a separable potential in §7]. Here 



justification means that non-stationary solutions of (1.1) can be approximated by CME dynamics over 
long but finite intervals. Given the analysis below, this result of [T3] can be immediately transfered to 
our non-separable case. J 



2 Band structure and Bloch functions 

Let uJn{k),n £ N, denote the spectral bands and Un{k;x) the corresponding Bloch functions of the 
operator L := —A + V{x), where k runs through the first Brillouin zone = (—1/2,1/2]^. This 
means that {ujn{k) , Un{k; x)) is an eigenpair of the quasiperiodic eigenvalue problem 

Lunik; x) = ujn{k)un{k] x), X e := [0, 27r)^, 

un{k- (27r, X2)) = (^^^^^Un{k; (0, xs)), tx„(fc; (xi, 2^)) = ^^^''^Un{k] (xi, 0)). 
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The Bloch functions Un{k;x) can be rewritten as 

Un{k;x) = e^^'^pnik; x), where Pn is 27r-periodic in both xi and X2, and fulfills (2-2) 
L{k;x)pn{k;x) := [{idx^-ki)'^ + {[dx2-k2)'^ + V{x)]pn{k;x) = u}n{k)pn{k; x) . (2.3) 

For each fc G the operator L(k; •) is elliptic and self adjoint in L^(P^), which immediately yields the 
existence of infinitely many real eigenvalues ujn{k), n G N with 0Jn{k) — )■ oo as n — t- oo. The spectrum 
of L equals UneN feeT^ '^n(^)! see Theorem 6.5.1 in ^14]. Moreover, via symmetries the ijJn{k) can be 
recovered from their values in the irreducible Brillouin zone Bq, see Fig. [2j From we also note 
that for a; G P-^ and n G we have 



Un{k; {xi + 2ni7r, X2 + 2n27r)) = e^™'^Un{k] x). (2.4) 




Figure 2: The first irreducible Brillouin zone Bq for the two-dimensional potential V . 

Gaps in the spectrum of L have to be confined by extrema of bands. Unlike in the case of the 
separable potential V{xi,X2) = Wi{xi) + W2{x2) the extrema of ojn within Bq do not have to occur 
only at k = T,X and M but may occur anywhere throughout Bq. Thus we need to solve (2.1) for all 
k G Bq. 



We choose the contrast rj of the periodic structure (1.3) so that two finite band gaps are open. Our 
computations show that this happens, for instance, at rj = 5.35, which we select. The band structure of 
L is computed in a 4th order centered finite-difference discretization. For reasons of tradition we plot 
in Fig. [3]the band structure along BBq. In Fig. |4]we plot the first few bands over B. Though not true 
in general [21], in our case the extrema of the first 6 bands fall on BBq. The dots in Figs.[3]and[4]label 
those band edge extrema which also mark gap edges. One of these extrema in Fig.js] (corresponding 
to 4 extrema in Fig. [4]) falls out of the vertex set {F, X, M}. We also label in Fig. [3] the first 7 bands 
cji, . . . , (xiy and the gap edges si, S2, • . • , S5. The edge values with six converged decimal places are 

si ^ 1.502064, S2 ~ 1.702299, S3 ^ 2.034433, S4 ~ 3.807113, and S5 ^ 3.832442. 



For any corresponding value of k each gap edge eigenvalue of (2.1) is simple because none of the 



edge-defining extrema belongs to more than one band. We now combine this with symmetries of the 
problem to find symmetries of the Bloch functions, which will be needed in the derivation of the CME. 
In the rest of this section we assume ||'Un(^; ■)IIl2(p2) = 1, where we are, of course, still free to multiply 
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Figure 3: (a) Band structure of L with r] = 5.35 along dB^. Red dots label band extrema at gap edges 
si, . . . , S5. (b) Detail in the second finite gap. (c) Detail near the edge S3 showing that 002 is not flat 
for k between X and M. 




Figure 4: Band structure of L with r/ = 5.35. On the right a detail of ^5 and we near the second finite 
gap- 



any mode Un by a phase factor e^", a G M. See also Remark 2.1 
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First, due to evenness of V{x) in both variables we have 



Un{{-ki,k2); (xi,X2)) 



'Un{{ki,k2); {2tT - Xi,X2)), 

'-Un{{ki,k2); (xi,27r - X2)), 



(2.5) 



Uni-ki, k2) = Unikl, -k2) = UJn{k). 



for some ai,a2 £ I^- Note that when {—ki,k2) = {ki,k2), where k = / reads "fc congruent to /" and 
means k = 1 + m for some m G Z^, a renormahzation of the phase cannot be used in general to obtain 
ai = because n„((fci, A;2); (tt, X2)) = Vx2 G P is possible. Similarly, when {ki^—k2) = {ki,k2), we 
cannot generally achieve 02 = because m„((/ci, ^2); (xi,7r)) = Vxi G P is possible. 
Next, the symmetry V{xi,X2) = V{x2,xi) implies 



Uni{kl,k2); (xi,X2)) = e'"'Uni{k2,ki)] (X2,xi)), UJn{ki,k2) =UJn{k2,ki). 



(2.6) 



for some a G M. Similarly to the case of symmetry (2.5), when ki = k2, one cannot, in general, apply 



renormahzation to achieve a = because Un{{ki, ki); (xi, xi)) = Vxi G P is possible. 



Finally, since L is real, Un{k]x) satisfies (2.1) with the factors in the boundary conditions replaced 
by e-^^'^^i and e''^'^'^^, xhus 



Un{-k;x) = Un{k;x), UJn{-k) = UJn{k). 



(2.7) 



Note that unlike in (2.5) and (2.6) no exponential factor appears in (2.7). This is because for the 
conjugation symmetry (2.7) such a factor e'" can be easily removed via multiplication by e~"*/^. 



Remark 2.1 If, e.g., (— fci, ^2) is not congruent to {ki, ^2), we can, for instance, multiply Un{{—ki, /C2); •) 
by e"^i and obtain Un((— ^1,^2); (xi,X2)) = Un{k] (27r — xi,X2)). However, one will generally not be 
able to simultaneously ensure also Un{{ki, k2)] {xi, X2)) = Un{{k2,ki); {x2,xi)) in (2.6) and therefore 
we stick to the factors in (2.5) and (2.6). J 



Let us consider implications of the above three symmetries (2.5), (2.6) and (2.7) for our exam- 



ple (1.3) and plot the gap edge Bloch functions in Fig. [sj Each edge si,S2 and S4 is attained only 
at a single extremum within B, namely at k = T,M and M respectively. The corresponding Bloch 
functions are tii((0, 0); x), iti((l/2, 1/2); x) and «5((l/2, 1/2); x) respectively, which are all real due to 
(|2.7[). The edge S3 is attained by extrema at A; = X and X' with the Bloch functions «2((l/2, 0); x) and 



U2{{0, 1/2); x). Referring to (2.6) only ii2((l/2,0); (xi,X2)) is plotted, which is again real due to ( |2.7[ ). 
Finally, the edge S5 is attained by 4 extrema, namely at A; = (kc, kc), {—kc, kc), {—kc, —kc) and {kc, —kc), 
where the numerically computed value, converged to 6 decimal places, is kc ~ 0.439028. The corre- 
sponding Bloch functions are UQ{{kc, kc); x), UQ{{—kc,kc);x), ue{{—kc,—kc);x) and UQ{{kc,—kc);x). 



Due to (2.5) and (2.7) and because kc ^ {0,1/2}, we can normalize the Bloch functions so that 
ue{{-kc,kc); {xi,X2)) = m{{kc,kc); {2Tr - xi,X2)), n6((/cc, -fee); (xi, X2)) = UG{{kc,kc);{xi,27r - X2)), 
ueii—kc, —kc)', (xi, X2)) = UQ{{kc, kc); {2-k — xi,2-k — X2)) = UQ{{kc, kc); (xi, X2)). Thus it suffices to plot 



only UQ{{kc, kc); (xi,X2)). In addition, (2.6) and the fact that UQ{{kc, kc); (xi,xi)) is not identically zero 
imply UQ{{kc, kc); (xi,X2)) = ue{{kc,kc); (x2,xi)). The Bloch waves ue{{kc, kc);x) and UQ{{-kc, -kc);x) 
are, therefore, symmetric about the diagonal xi = X2. 
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s,: u,((0,0);x) 



Sj: LJ,((1/2,1/2);x) 




Figure 5: Bloch functions at gap edges si, S2 • • • , S5. 



Fig. [5] shows that all the Bloch functions at si, S2, • • • , S4 are either even or odd in each variable. 



This actually follows from (2.5) and the fact that these gap edges occur at k £ T, = {T, X, X' , M}. As 
each coordinate of any A; G S is either or 2, the eigenvalue problem (2.1) is real and we can choose 



(11,0-2 £ {0,7r} in (2.5). The choice ai 
(2.5). Taking, for instance, ki 



02 = is, however, in general impossible as explained after 



|, we have 



ti„((l/2, k2);x) = ±u„((-l/2, k2); (27r - xi,X2)) = ±'u„((l/2, k2); (27r - xi, X2)) 

= ±e^27r|^^((^/2, k2); X2)) = TWn((l/2, ^2); {-Xl,X2)), 

where the second equality follows from 1— periodicity of Un in each A;— coordinate and the third equal- 



ity from the quasi-periodic boundary conditions in (2.1). Similarly, we get Un{{0,k2);{xi,X2)) = 
±ii„((0, /C2); (—2:1, X2)). Therefore, we have the following 

Lemma 2.2 Suppose V{x) is even in the variable xj for some j E {1,2}. If kj G {0, 1/2} and iOn{k), 
as an eigenvalue of (2.1) has geometric multiplicity 1, then Un{k;x) is either even or odd in Xj. 
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3 Formal asymptotic derivation of Coupled Mode Equations 



Gap solitons in the vicinity of a given band edge are expected to be approximated by the Bloch waves 
at the band edge modulated by slowly varying spatially localized envelopes. The governing equations 
for the envelopes, called Coupled Mode Equations (CMEs), can be derived by a formal asymptotic 
procedure. Here we are interested in gap solitons E{x,t) = </>(x)e~"^* with uj = uj^ + e^Jl, < e <C 1, 
where 0;=,, is an edge of a given band gap of a fixed (0(1)) width, and has a sign chosen so that uj lies 
inside the gap. The leading order term in the asymptotic expansion of the spatial profile (p is expected 
to be 

N 

<P{x) ^e^Aj{ex)un^{k^^^^;x), (3.1) 
i=i 

where {un {k^-^''',x)\^ are the Bloch waves at a; = and (2.4) is used for x P^. We assume: 
Assumption A.l The band structure defined by (2.1 ) has a gap with an edge (lower/upper) defined 
by < < oo extrema (maxima/minima) of the bands a;„ ( /c) . The extrema occur for bands {k) , j = 
1, . . .N at the corresponding points k^^^ G B. 

Assumption A. 2 The quadratic form dl_^u}nj{k^~'^)x'^ + 2dk^dk2ijJnj{k^-'^)xy + dl^u}nj{k^~'^)y'^ defined by 
the Hessian of a;„. at A; = k^^') is (positive or negative) definite. 



Remark 3.1 Assumptions A.l and A. 2 are satisfied by the potential (1.3) with rj = 5.35 at all the 



gap edges si,...,S5. J 

Remark 3.2 The definiteness in A. 2 ensures that the extremum of ajn^ at k = k^^^ is quadratic and 
that the resulting CMEs are of second order. Unlike in the separable case [19] it is possible that 
dkidk2i^njik^''^) 7^ 0, which then leads to CMEs with mixed second order derivatives. J 

Remark 3.3 The Bloch waves Unj{k^-'^; •),j=l, . . . ,N defined by the extrema are called "resonant".] 



Remark 3.4 The approximation (3.1 ) with the same e-scaling applies also to gap solitons in an ©(e^)- 
wide gap which closes at = cj,,, as e — >■ in such a way that the plane w = a;^, at e = is not intersected 
by any band but is tangent to bands at extremal points. Um , ■ ■ ■ , Unj^ are then the resonant Bloch 
waves at w = a;=K at e = 0. Such a case was studied in [13j for a separable periodic potential. 

The above discussion is not limited to the case of the Gross-Pitaevskii equation but applies to 
general differential equations with periodic coefficients, as it depends only on the band structure. A 
typical example is Maxwell's equations with spatially periodic coefficients. J 



We now give the derivation of CMEs under the assumptions A.l and A. 2. For the example (1.3) 
with 7] = 5.35 we first review the derivation in physical variables (p{x) near u = S3, then comment on 
an obstacle for this calculus near to = s^, and therefore present a derivation in the general case in the 
so called Bloch variables which avoids this obstacle. Finally, we apply this general procedure to all the 



five gap edges of the example (1.3). 
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3.1 CME derivation in Physical Variables 

The ansatz in physical variables is 



AT 



e,/>(o)(x) = e^^j(y)u„^.(/e(-^');x), u = + e'^Q, y = £x, < e <C 1. 



(3.2) 



To review the derivation of the CMEs we choose = S3 for the example (1.3) with rj = 5.35 



3.1.1 CMEs near the gap edge = S3 

At the edge S3 we have N = 2,ni = n2 = 2,k^^^ = X and k^'^^ = X' , i.e. the two resonant Bloch 



waves are vi{x) := U2{X;x) and ^2(2;) := U2{X';x). Using (2.4), Lemma 2.2 and (2.7), we have that 
vi is odd and 27r— antiperiodic in xi and even and 27r— periodic in X2- Opposite symmetries hold for 



V2- Moreover, (2.7) implies that vi and V2 are real. We normalize the Bloch functions vi^2 over their 
common period [— 27r,27r]^ so that ||fj||L2Q_27r,27r]2) = 1, j = 1,2. 



Substituting (3.2) in (1.2) leads to a hierarchy of problems at distinct powers of e, each of which 



we try to solve within the space of functions 47r-periodic in both xi and X2, invoking the Fredholm 
alternative (see e.g. chapter 3.4 of |35| ) where necessary. At 0{e) we have the linear eigenvalue problem 
[L - S3]vj{x) =0, j = 1, 2. At ©(e^) we have 

[L - S3](l)^^^ = 2 {dy^Aidx^Vl + dy^A2dxiV2 + dy^AAiVl + 5y2^29a;2?^2) • 



By differentiating the eigenvalue problem (2.1) with respect to kj,j € {1,2} and evaluating at n 
2,k = X = (1/2,0), we find that 



[L - S3]v[^'\x) = 2dx^vi, 



(3.3) 



and similarly [L — S3]v2 ^ (x) = 2dx V2, where 



v^^^\x) = -i(9..p2(X;x))e'^-" and v^'^'^x) = -i(9. p2(X'; x))e 



JX'-x 



are called generalized Bloch functions [29]. Thus 4>^^^ = dy-^Aivl^^^ +dyj^A2V2^^^ +dy^Aiv[^^^ ^dy^A 



,(^1) 



,(^2) 



(x ■) 

(3.3) implies that Vn ' {x) is odd/even in Xj if Vn{x) is even/odd in Xj respectively. 
At ©(e^) we obtain the CMEs. We have 

[L - S3](/)(2) = 9.{AlVl + A2V2) + /^y^^y^AlVl + l\y^^y^A2V2 

+2 [dlAA.vt''^ + dlA2dx,v^r^ + d^AAA"^ + dlA2dx,v^P^ 

^dy^dy^Axdx^v'x^'' + dy^dy^A2dxiV^2^^ + dy^dy^Aidx^v^i''^ + dy^dy^A2dx2V^2 
Ei=i l^iP^i^^j + 2|Ai|2A2wfu2 + 2\A2\'^ Alvlvi + AlA2vlv2 + AlAivlvi , 

and the Fredholm alternative requires the right hand side to be L^{—2-k, 27r] ^-orthogonal to vi and V2, 



2V2 



(^2) 
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the two generators of Ker(L* — S3). Taking the inner product, we see that the terms 1, ^2) and (f2, f 1) 
in the inner product vanish due to orthogonahty of Bloch waves. Many additional terms vanish due to 
odd or 27r-antiperiodic integrands (in at least one variable). Namely, in the inner product of the right 
hand side with vi the integrals {dxiV2^^\ vi) , {dx2V2^'^\ vi) , {dxiv[^^\ vi) and {dx2v[^^\vi) vanish due 
to odd integrands and the integrals (t>2,wi), {viV2,vi), {viV2,vi) {dx-^v^^\vi) and {dx2V2^^\vi) due to 
27r— antiperiodic integrands. An analogous discussion applies for the orthogonality with the respect to 
V2- The remaining terms have to be set to zero, which leads to the CMEs for the envelopes Ai and A2: 



nAi + aid^^Ai + a2dlAi - a [tiI^iI^^i + -f2i2\A2\^Ai + AjAi)] =0, 
nA2 + a2dl^A2 + aidl^A2 - a + -i2{2\Axf A2 + A\A2)\ =0, 

I vidxj^v^i dx, a2 = 1 + 2 / / vidx2Vi dx, 

27r J -2iT J-2tt J-2n 



(3.4) 



27r f2n f2n p2TV 

I vfdx and 12=1 / 



3.1.2 CMEs near the gap edge 55 

At = S5 we have = 4. The resonant Bloch waves are vi := UQ{{kc, kc); x), V2 := UQ{{—kc, kc); x), 
^3 := UQ{{—kc,—kc);x) and V4 := UQ{{kc, —kc); x). Analogously to ^3.1.1 the asymptotic expansion 



needs to be carried out in the space of functions periodic over the common period of . . . , ^4. The 
Bloch functions are then pairwise orthogonal over this domain. However, if kc is not rational then 
the Bloch waves are not periodic but only quasi-periodic. Therefore, unlike in the case of a separable 
V{x) [13], where always kc £ {0, 1/2}, in the non-separable case in general the derivation in physical 
variables is impossible. 

3.2 CME Derivation in Bloch Variables (p(k; x) 



An alternative to the derivation in ^3.1 is to transform the problem to Bloch variables. The advantage 
is that the linear eigenfunctions are then all 27r— periodic in each x— coordinate. The orthogonalization 
domain is, therefore, always P^. 

3.2.1 General Case 

The Bloch transform T is formally defined by 

4>ik; x) = iT<P)ik; x) = ^ e'"'--4>ik + m), 0(x) = (T^'m^) = e''-^4>{k; x)dk, (3.5) 



where 4'{k) denotes the Fourier transform. T is an isomorphism from i7*(M^, C) to L^(T^, H^'^"^ 



2 



L2(T2,/i-«{p2,c)) 



j^2 \\(t>{k'-, ')\\'\{s^^2^dk, cf., e.g.,[3T], and by construction we have 

(j){k; (xi + 27r, X2)) = <^(A;; (xi, X2 + 27r)) = 0(A;; x), (3.6) 
0((A;i + l,fe2);x) = e-'^'i(^(A:;x), 0((A:i, ^2 + 1); x) = e-'^'2<^(fe; x). (3.7) 
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Multiplication in physical space corresponds to convolution in Bloch space, i.e., 



J2 



4>{k — I; x)ijj{l; x)dl =: [4>*b 2;), 



(3.8) 



where (3.7) is used \ik — I ^ T^. However, if g is 27r— periodic in both xi and X2, then {T{gu)){k; 



g{x){Tu){k;x). 

In order to choose a suitable asymptotic ansatz for (p{k;x), note first that the Bloch transform T 
of the ansatz (3.2) for e(j)^^\x) is 



:4>^'^\k; 



1 ^ / 

i=i mez2 \ 



k — k^^^ + m 



with k £ T'^,x £ F'^. As is localized near p = 0, we approximate Aj 



(3.9) 



by XD, (^ + 



m)Aj 



and 



. / k-k'-^^+m 



where xDj (k) is the characteristic function of the set 



1 2 
- < r < -. 
3 3 



(3.10) 



(3.11) 



The reason for (3.11) will be explained in ^4.2 



Below we will also use periodically wrapped versions Dj of these neighborhoods, i.e. 

Dj := {fc G : \k-k^^^\ < e' modulo =} 



(3.12) 



where 'modulo =' means equal modulo 1 in each component, see Fig. [6] for an example. 

Note that k + m £ Dj with A; G is possible only for m G {m G : —1 < mi, 7712 < 1}. We 
define the set of m— values for which k + m G Dj for some A: G by Mj := {m a 1? : k + m £ 
Dj for some A; G T^}. In fact, for small e only the following cases occur: Mj = {( ) , ( 

)} if k^^ = 1/2 

and fcO) / (1/2, 1/2), = {( g ),(?)} if k'i^ = 1/2 and fcO) / (1/2, 1/2), M, = {( g ),( J ),(?),({ )} 
if k^i) = (1/2, 1/2), and Mj = {(§)} if k^^) G int(T2). 

Thus we are lead to the following asymptotic ansatz in Bloch variables 



^{k-x) = -^'^^^\k;x)+i^^^\k-x) + ei)^'^\k-x) + 0{e^), 



i^^^\k;x) 



N / 

^p„^.(A:W;x) XD,{k + m)Aj\ 
j=i jTteM, \ 



k - k^i^ + 



m 



(3.13) 



< e < 1. 



The periodic part of the Bloch functions is normalized so that \\pnj{k', •)||/^2(p2) = 1. 



12 




X' 




X ki 



Figure 6: Sets Dj and Dj for k^^'^ = X and A;^-'^ = X' (as in the example ( 1.3 ) with r] = 5.35 at co^, = S3). 



The difference between the leading order terms in (3.9) and in (3.13) is 

N 

e 



1 ^ 



X 



(3.14) 



with 



hj{k;x) = ^pnj{k^^^;x) 



{I - XD.ik + m)) A 



3 \ e 1^ 



me Mi 



k—k^^^+m \ im-x 



(3.15) 



We now estimate ||/ij||^2C]['2^jys(p2)). In the first sum in (3.15) we have \k + m — k^^^\ > while in 



the second sum \k + m — fe^-'^l > 1 because k + m ^ Dj for all A; G if m G \ Mj. By the triangle 
inequality and the substitution p = {k — k^^^ + m)/e we obtain 



\h 



'jllL2(T2,_ff^(p2)) 



< 



E IK(A;(^'V)e^"^-|li.(p.) / 



pe(T2-fc(j)+m)/e 



\Aj{p)\'^dp 



+ E IK(^^'V)e^'"-| 



|p|>£'- 



\Aj{p)\^dp + 



pe(T2-fcO)+m)/e 
i,(p)|2dp , 



where the regularity of Pujik^^'^; ■) is guaranteed if F G H^^^{^^). By rewriting the right hand side 



as C 



out of the integrals, we have 

ll^illL2(T2,//=(p2)) < C(e''(^~'') + e'')||ij||i2(R2) < Ce'*(^"'')||lj||i2 



and taking the supremum of ( 1 + |p| ) 

(3.16) 
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For r < 1 we thus have that e '^il)^'^\x) approximates £4>^^\x) up to 0{£^^^ ''^) in the H^{^) norm. 
Because ||jys(]R2) = Oil), this approximation is satisfactory. 



Applying next T to (1.2) yields 



L — uj 



• + a (j)*B 4> *B 4> = 



L{k;x) = L{k^^^ 

= L{k^-^'^ -m; x) - 2e 



m + ep^^'^'> ; x) 



(i) 



Pi +P2 



(3.17) 



on {k;x) G x P^, where we recaU from ([2^ that L{k;x) = (id^^-kif + {id^^-k2f + V{x). 
Setting := ^+"7^''' , we have 



(3.18) 



Substituting (3.13) in (3.17) and using (3.18), we obtain a hierarchy of equations on a; G P^,A; G 
such that k + m G Dj,j G { 1 , . . . , . Note that the combination of /c G and k + m ^ Dj impHes 
m G Mj. The following hierarchy is thus for each (j, m) G {1, . . . , A^} x Mj. 



©(e-i) : ij(p(^'™)) Z(A:(J') -m;x) -w, (p„. (A;(^'); x)e'™-^) = 0, 



which is equivalent to Aj{p^^''^^)e^"^'^ L{k^^');x) — oj^ 

t^* = UJnj{k^^^). 

0(1) : Z(A;(^) — m;x) — 



PnAk^^^;x) = and thus holds by definition of 



2i,(p(^'")) p(^''")(i9,,-A;l^'^ +mi) +p^^''"^(i5,,-4^'^ +m2) {pn,{k'^'^;x)e''^-^) 



2ij(p(j''"))e^™-^- 



pS^''")(ia.,-A:i^-))+p(^-'")(ia.,-4^)) 



,(fc(^-);x) 



for fc+m G -Dj, j G {1, . . . , N}. To solve this, we note that by differentiating [L{k; x)—ujnj {k)]pn^ {k; x) 
with respect to ki,l G {1, 2} and evaluating at k = k^^^ — m, we obtain 



L{k^^^ -m;x) - dkiPnj {k^-''' -m;x) = 2(i5^,-A;p^ + mi)pnj{k^^^ - m; x) 



(3.19) 



Since Pn{k — m; x) = e^"^'^pn{k; x) due to (3.7), we get for A; + m G Dj 



^^^\k-x) = J]pp'™U,(p(^'™))e'™-^9fc,p„^(fc(^');x). 



(3.20) 



1=1 
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0{e) : We have 

; X 



2 



(3.21) 



using V^^i) from (|3^. 

The nonlinear term has the form 



Gj{k; x) ■.=-^XD, {k + *b 4^^'^ *b P^){k; x) = -^xd, {k + m) 



N 



*B *s Co 



.0 = 1 



TV 



N 



N 



+2 ^ *B C/3 *BS,a+ X] *s *s + X] (.a *B (./S *B 



a,/3=l 



where = ^a{k;x) := PnAk^'^'^-x)Y,m&M^XD^{k + m)Ac 



a,/3,7=l 



A;+m— \ „im- 



(3.22) 



fc— m+fct") \ im-x 



The last sum or the three last sums in (3.22) 



are absent if = 2 or = 1 respectively, ^q, *b *b S,^ consists of terms of the type 



5no,(A:;x) =ei("+°-«)-VJA;(-);x)p„^(fe('^);x)p^(feW;x) / / XD„(A;-r+n)A 



e 



X XDAr-s+o)Ai3 



-g+o-fc('^) 



T2 T2 

x-D^(s-g)^7 



(3.23) 



ds dr 



with n G Mq, o S Mg and g G M^. Clearly, the integration domains can be reduced to r € D2e'-{k^^^ - 
— o + q) and s G Df.r[—k^'^'^ + g). The changes of variables s := (s + k^'^^ — q)/e, and f := 
(r - /c^/^) + /c('>) + o - g)/e yield 



9noq{k] x) — e e '- ^ Vndk^ ] x)pni^{k^^'^ ] x)pn^[k^'^^ ; x)x 

J ( fc-(fc(")+fc(^)-A:(^))+n+o- g 



XD, 



r X 



2£' e e' £ 



(3.24) 



i« ( ^-^^'"'+''''7'^^'^+"+°-" - r ) XD (r - 5)A/3(r - .^XD {s)A,{s) d5df , 
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where D^r-i = {p£M? : \p\ < e''^^}. 

Only those combinations of (n, o, q) which produce nonzero values of all the three characteristic 
functions in (3.23) for some k,r,s S are of relevance. Due to X-D^is—q) we, therefore, require 
g-fcW G T2 = [-l/2,l/2]2, which ensures that s — q £ —D^ is satisfied by some s S for any e > 0. 
The first condition is, thus. 



so —q-k^^) ET2. 
Due to xopi'i'—s+o) we get the condition sq — o + k^^^ G T^, i.e., 

ro := SQ-o + k^P'^ E P. 

Finally, XDa(^~'''+^) enforces ro — n + /c^") G T^, i.e., 

feo := ro - n + A:^") G T2. 



Statements (3.25), (3.26), and (3.27) form the necessary condition 

sq:= q- k'^'^'^ G T^, ro := sq 
for (3.23) (and thus ( 3.24[ )) not to vanish. 



o + A;(^)gT2, and ko :-- 



ro — n 



(3.25) 



(3.26) 



(3.27) 



(3.28) 



Another condition on {n,o,q) appears due to the factor XDj{k + fn) in Gj. From (3.24) it is clear 
that Qnoq is Supported on A; G Df^r (A;(°) + k^l^^ — k^^^ — n — o + q). The factor xDj {k + m) thus annihilates 
all terms gnoq except those for which 



— n — o + q 



m. 



(3.29) 



If (3.29) is satisfied, (3.24) becomes 

9noq{t,x) = e^e'("+°-'?)-V.(A;^"^;^)Pn,(A:^''^x)]5^(A;(^);x)x 



/ 



XD I ^=^5^ - r I X 



-D„r-in 



A:-fc(j)+m 



T2+fc(7)-q 



(3.30) 



As a result, the term AaApA^ will enter the j— th equation of the coupled mode system provided 
there exist n G Ma,o G and g G such that (3.28) holds and such that (3.29) holds for some 
m G Mj. Let us denote the set of {n,o,q) that satisfy (3.28) and (3.29) by Aa,p,'y,j,m- 

The sum of the terms (3.30) over (n, o, q) G Aa,i3,'y,j,m yields a double convolution integral over the 
full discs f G D2i;r-i and s G D^r-i, i.e.. 



(Ca *B *B e^)(A:; x) = eV('='"'+'=''''-'=*^'-'='^'+-)-V. (A:^"); (A:(/^); x)^^ 



'Do 



-1 ^i^.r-l 



r\Ar_ 



k—k^^^+m 



n XD^r-i{r - s)Ai3{f - s)A^{s) dsdf, 



(3.31) 
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where e'("+°-5)-^ was replaced by ei(fe(")+fc(«-fc(^'-fc(^>+m)-x 

We return now to equation (3.21 ) for ^^^^ on A; G (Dj — m) PlT^. Its solvability condition is L-^(P^)- 
orthogonality to Ker(Z(/c(-') — m;x) — w*) = span{U;Pri, (/c^-'); a;)e^™'^ s.t. a;„,(A;(-^^) = a;*}. Clearly, the 
dimension of the kernel is at most N. The value is attained if k^^^ = . . . = k^^\ 

In the linear terms in (3.21) the factor e'"*'^' is canceled in the inner product with pn^{k^^^ ; x)e^'^'^ 
so that the same solvability condition holds for all ni. The range of is a different section of the 

disc D^r-i for each m. The section is an (l/|Mj|)-th of the full disc so that these \Mj\ conditions build 
one equation in p S D^r-i. 

The resulting N equations are CMEs in Fourier variables p G D^r-i: 

Oi, - {^-d^un, ik^'^)pi + Idl^n, ik^'^)pl + dk.dk.u^n, ik^'^)piP2^ A, - = 0, (3.32) 

i e {1, . . . , iV}, where AA^ (p^^''") ) = (Gj(ep(j'™) + k^^^ - m; ■),Pn0^^; •)e""-)L2(p2). 

For sufficiently smooth Aj we can neglect the contribution to Aj from p G \ D^r-i or, for 

simplicity, assume that the Aj satisfy (3.32) also there. Equation (3.32) is then posed on p G M^. 
Performing the inverse Fourier transform yields the CMEs 



nAj + 



0. 



(3.33) 



The structure and coefficients in Mj for our example (1.3) will be discussed in ^3.2.2 



In order to make the discussion of the asymptotic hierarchy complete, we need to mention the part 
of the /c— domain outside the neighborhoods of k^^\ For A; G such that k + m G T'^ \ Dj for all 
m G Mj we have [L(A;(J') -m;x) - w*] ^^("^(A;; x) = forn = 1,2 so that Tp^°\k; •) = ip^'^^k; ■) = for 
such k. 



The appearance of second derivatives of the bands conj in (3.32) is due to the following 
Lemma 3.5 For any l,m G {1, 2} 

dkAn.^n,{k^'^) = 25i^ - 2{{id,^-kii^)dk,Pn,{k^^^; •) + {id,-kl^^)dk^pn^{k^^^; •),Pn,(fc(^'); •)>L2(P2), 

where 5im is the Kronecker delta. 

Proof. This follows from differentiation of [L{k; x) — a;„^. (A;)]p„^. (/c; x) = w.r.t. k. □ 
As the next lemma shows, for even potentials V{x) the mixed derivatives of ujnj are zero whenever 
u}n-{k^^^) has geometric multiplicity one and the extremal point A;^-'^ coincides with one of the vertices 
of the first irreducible Brillouin zone or of its reflection. 

Lemma 3.6 Suppose V{x) is even in xi as well as in X2. Then dk-^dk2'^nj{k^^^) = dk^dk-^ojnjik^^'^) = 
if k^^^ G S = {r,Ar, X',M} and provided ojn\k^^^) has geometric multiplicity 1 as an eigenvalue of 
(23. 

Proof. Take l,m G {1, 2}, I ^ m. As idx^—km is self-adjoint, we have 

{{\dx^-k'i)dk,PnM'^; ■),Pn,{k^'^; •))l2(P2) = {dk,Pn0''^-r), {id,^-4i>)pn,{k'^^^ ; ■)) L^^2y (3.34) 
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Based on (2.2) we have iidx^—km)pnAk^-'^',x) 



-le 



"'^'^■''d,^UnAk(^^;x). Next, dk,PnAk'^'^;x) 



le 



ifc(j'-a;^{^^i)|'^(j). ^-j^ where Vnj\k^^^;x) is the generahzed Bloch function [29] solving 



[L - u:^]u = 2dxiUnAk''^^]x) 



u(27r, X2) 



(7) 

j^^"'^! u(0,X2), n(xi,27r) 



tx(xi,0), (3.35) 



analogously to (3.3). The inner product in ( |3.34 ) thus becomes {—Vn-\k'^^'>\-),dx,^Unj{k^-'^;-)) 



/L2(p2). 



Because A;^-'^ G S, Un^ is even or odd in xi (Lemma 2.2). From (3.35) it is clear that Vn^ {k^^^; x) has the 
opposite symmetry (odd or even respectively) in x;. Thus, the integrand is odd in xi and the integral 
vanishes upon shifting the integration domain to [— tt, tt]^. □ 



3.2.2 CMEs for the Example (1.3) 



We now calculate the explicit form of the CMEs (3.33) in the vicinity of the five gap edges in the 



example (1.3) with 77 = 5.35. It turns out that only few terms are nonzero in the nonlinearity Nj for 
this case. Of special importance is the edge = S5, where k^^"* ^ S and, indeed, dk^dk20Jnj{k^^'') / 0. 

the functions 



In order to numerically evaluate the coefficients dkidk^0Jnj{k^^^) given in Lemma 



3.5 



dkiPnj{k^^^',x) have to be computed. They are solutions of the singular system (3.19) but as the right- 
hand side is orthogonal to the kernel of L{k^^^; x) —oj^,, the BiCG algorithm can be used as long as the 
initial guess is orthogonal to the kernel. We work in a 4th order finite difference discretization and use 
an incomplete LU preconditioning for BiCG. 



CMEs near = sii Only one extremum defines the edge w* = si, namely the minimum of the band 
ui at k = r. Therefore, N = l,ni = I and k^^'> = F. Because fc^^) G int(T2), we get Mi = {(0,0)^}. 
Thus 



[n + a{dl+dl)]A-aj\A\'A = 0, 

where a = ^5|^^i(r) = |5|^c^i(r) and 7 = (pi(r; .)2,pi(r; •)2)^2(p2) = \\pi(T;-) 
in a holds due to (2.6). The numerically obtained values are a w 0.62272 and 7 



IL4(P2)- 

a 0.048029. 



(3.36) 
The identity 



CMEs near uj^ = S2'- Here the linear problem is characterized by = l,ni = 1 and k^^^ = M = 
(1/2, 1/2) and we get Mi = {(0,0)^, (1,0)"^, (0, 1)"^, (1, 1)^}. The resulting CMEs have the form 
We determine next the coefficient of the nonlinearity In (3.24) we have a = /3 = 7 = 1 and 

^(") + k^P) _ - fcO) = (0,0)^. We carry out a straightforward sweep through all the possible 
combinations {n,o,q,in) (performed using a Matlab script) to determine those that satisfy (3.28) and 



(3.29). As a result we have 



™.M.(:.1)M§): (..,°..)^^{(H)-(l8)-(8l)-(!!)} 
».AM:.2)MJ): »• ^ { (g §) > (I ») ■ (S l) ' (I 1 ) } 
™.AM:.3)M?): »• ,f ^ { ( ^ ) • ( 1 ») ■ (g 1) ' ( ! I ) } 
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m 



AM:,4) = (!): („, „, ,f e { ( J J ) , ( | J) , ( J | ) , (| | ) } 



where we have used Mj{:, I) to denote the /— th vector in Mj. 

The CME coefficients are thus a = \dl^uji{M) = \dl,ui{M) and 7 = (pi(M; •)^)l2(p2) = 

\\pi{M; •)lli4(p2)- The identity in a holds due to Numerically, a w -1.971217 and 7 w 0.076442. 



2, = X and fe^^) = x'. We have thus Mi 



CMEs near uo^ = S3: Here N = 2,ni = n2 
{(0, 0)^, (1, on and M2 = {(0, Of , (0, 1)^}. 

For Mj we sweep again through all the possible combinations (n, o, g, m) for both j 
The results are summarized in Table [TJ 



1 and j = 2. 



term 
in Afi 



(n,o,qy satisfying (3.281 and (3.29 1 



m = Mj{:, 1) 



M,(:,2) 



coefficient of 
the term in aNi 



\A2\^A2 



TIT 



(!]) 



1/2 
-1/2 
-1/2 

1/2 



-1/2 
1/2 





1/2 
-1/2, 



(p2(X,.f,P2(X,.)^) 



(p2(X',.f,p2(X',.)^} 





2{\p2{X,-)\\,\p2{X',-)\') 

2(b2(X,.)|Mp2(X',.)H 







AM2 



1/2 
-1/2, 









(e'(^■-^'"■p2(x,.)^P2(x^.f) 



-1/2 
1/2 



(e'(-l■^'^■>2(x^.)^P2(x,.)^) 



Table 1: Calculation of the nonlinearity terms for the CME near = S3. 



The resulting CMEs are 

[f] + aidl^ + a2dl] Ai-a [ji\Ai\''Ai + 72(2|^2|'^i + AjAi)] =0, 
[n + a2dl^ + ai^^J A2-a [71IA2PA2 + 72(2|^i|2A2 + ^?^)] =0, 

where 

ai = ^a2^u;2(X) = \dluj2{X'), 02 = ^9|^a;2(X) = \dluj2{X'), 

71 = (i'2(-'^; •)^^'2(-'^; •)^)l2(p2) = (P2(-'^'; •)^^'2(-'^'; •)^)l2(p2) 

= l|p2(X;-)|li4(p2) = ||p2(X';.)|li4(p2), 
^2 = (e^(i'-i)"p2(^;-)',P2(^V)'>L2(p2) = (e'(-i'i)"p2(XV)^P2(X;-)=^)L2(p2) 
= (|P2(X;-)|Mp2(^V)P>l2(p2). 
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The identities in ai, 02 and 71 hold due to (2.6). The equalities in 72 yield 72 € M and follow from the 
fact that U2{X,x) = e^^^/'^p2{X;x) and U2{X',x) = e^^^^'^p2{X';x) are real. In detail 



/ e''^'p2{X;xfe-'''^P2{X';xfdx= [ U2iX; xfu2{X'; xf dx 

Jp2 JpZ 



U2{X;xyu2{X';xYdx= / n2(X'; x)^n2(X; x) dx = / e'""^ p2{X' ; xy e-'""' p2{X ; x) dx 



The CMEs (3.37) are thus identical to those in Sec. 3.1.1 derived in physical variables. Numerically, 
ai w 2.599391, a2 ~ 0.040561, 71 w 0.090082, and 72 w 0.003032. 

CMEs near a;* = 54: Here = l,ni = 5 and k^^^ = M. This case is completely analo- 
gous to 0;=, = S2. The CMEs are ( |3.36D with a = |9|^W5(M) = ^dlur,{M) ^ -0.300655 and 
7 = (p5(M; ■f,P5{M; •)2)i2(p2) = ||p5(M; •)||i4(p2) ~ 0.039755. 



CMEs near = 55: Here iV = 4, ni = n2 = ^3 = 714 = 6, = {kc, kc), k^'^'^ = {-kc,kc),k''^^ = 
{-kc,-kc) and k^^^ = {kc,-kc), where kc ^ 0.439028. This is an important case in our example 
because k^^'^ ^ S here. Note that because k^^^ G int(T^) for all j € {1, . . . , N}, we have Mi = . . . = 
M4 = {(0,0f}. 



We start with the last two sums of G (see (3.23)). Terms of the type *b 6 *b Cm (the third sum in 
G) do not contribute to the CMEs because 2k^''^ — k^"^^ is not congruent to any k^^\j G {1, ... ,4} for 
any choice of m G {1, . . . , 4}, / 7^ m. For example, 2k^^^ — k^"^^ = (3A;c) k^), which is not congruent to 
any k^^^ since kc ^ {0, 1/2}. Only four terms of the type *b im *b in (the last sum in G) contribute 
to the CMEs, namely ^2 *s ^4 *s to the equation for k G Di, ^3 *b C4 to the equation for k £ D2, 
(,2 *B '^4 *B to the equation for k & D3 and ^1 *b ?3 *b ^2 to the equation for k G D4. This is because 

fc(2) + fc(4) _ ^(3) ^ kW^kW + ^(3) _ ^(4) ^ fc(2)^ ^(2) ^ ^(4) _ ^(1) ^ ^(3) ^(1) ^ ^(3) _ ^(2) ^ ^(4)_ 

The other terms in the last sum in G do not contribute. As an example, k^^^ + k^"^^ — k^^'^ = {kc, 'ike)- 

Another consequence of k^^^ ^ S is that Lemma 
may appear. The system of CMEs thus becomes 



3.6 



does not apply and mixed derivatives of Aj 



= [n + ai{d^^+d^y^) + a2dy,dy,] Ai 

- a 7i|^i|^^i + 2(72 (|^2p + \A^?) M + 7i|^3p^i + 72^2^4^) 
= [J] + ai(9^,+52j - a2dy,dy^\ A2 

- a 7l|^2|^^2 + 2(72 (|Aip + 1^31^) A2 + 7l|^4p^2 + ^2^1^3^4) 

= [r? + oixipl^^dl^) + oc2dy,dy^ Aj, 

- a 7l|^3p^3 + 2(72 (|^2P + |^4P) ^3 + 7l|^lP^3 + 72^2^4^) 

= [J7 + ai(9',+52j - a2dy,dy^\ ^4 

7l|^4|^^4 + 2(72 (|Ai|2 + |^3|2) A4 + 7l|^2p^4 + ^2Al^3^) 



(3.38) 
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where 



a2 

71 

72 

71 

72 



\d\i^^(k,{{-\T \-\Y)) = \dl^^^(k,{{-\y> for any m,n,p,q G {0,1}, 
dkidki^eikc, kc) = dk^dk2(^Q{-kc, -kc) = -dk^dk^i^ei-kc, kc) = -dk^dk^ueikc, -kc), 

{\pe{k,{{-ir, i-ir); ■)\^ ipeHk^a-ir, {-in, ^^^^(p.) = \\pem, k,)- 

for any m,n £ {0, 1}, 

{\P6{{-kc, kc); OP, \P6{{kc, kc); •)P>L2(p2) = {\P6{{kc, -kc); OP, \Pe{{kc, kc); •)P>L2{p2) 
{\P6{{-kc, -kc); OP, \pe{{-kc, kc); •)P)l2(p2) = {\pe{{kc, -kc); OP, l?'6((-fcc, -A^c); •)P)l2(p2), 
(|P6((-A:c, -A;c); OP, |P6((A;c, kc); •)P)l2(p2) = {\pe{{kc, -kc); OP, Ipeii-kc, kc); •)P)l2(p2), 



{Peii-kc, kc); ■)p&{{kc, -kc); ■),p%{{-kc, -kc); ■)p&{{kc, kc); ■)) 



The identities in a i, 02 and 71 are due to (2.5) and the identities in 72 due to (2.5) and (2.6). 

Moreover, 71 = 71 and 72 = 72 due to (2.7). This also imphes 72 = 72- Using these identities, we 
arrive at the system 

+ ai{dl^+dl^) + a2dy,dy,] Ai-a [71 (|Aip+2|A3p) Ai + 272 ((|^2p+| A4p)^i + ^2^4^)] =0, 
+ ai{dl^+dl^) - a2dy,dy,] A2-a [71 (| ^2^+21 A4P) A2 + 272 {{\Ai\''+\A^\'')A2 + ^1^3^)] =0, 
[17 + ai{dl^+dl;) + a2dy,dy,] A^-a [71 (|A3p+2|^i|2) ^3 + 272 ((|^2p+|A4p)^3 + ^2^4^)] =0, 
[Q + ai{dl^+dl;) - a2dy,dy,] A^-a [71 (|A4p+2|^2H ^4 + 272 ((|^iP+|A3p)^4 + ^1^3^)] =0. 

(3.39) 



The numerical values of the coefficients are ai 
72 « 0.029926. 



6.051248,02 ~ 0.096394,71 w 0.039118 and 



4 Justification of the Coupled Mode Equations 

If = a;=K + e^U, is in the band gap, then families of solitons, i.e., of smooth exponentially localized 
solitary wave solutions, are known for many classes of CMEs [33]. However, as already noted in the 
introduction, the formal derivation of the CMEs in ^ discarding some error at higher order in e. 



does not imply that localized solutions of the CMEs yield gap solitons of (1.2). For this we need to 
estimate the error in some function space and show persistence of the CME solitons under perturbation 
of the CME including the error. We proceed similarly to [ISJ. However, as function space we choose 
H^{^) with s > 1, in contrast to F^^L].(^) in [13]. The latter is possible in the separable case 
but there is the technical obstacle of the extension of |13l (3.7)] to the nonseparable case. On one 
hand, L^(M^) in Fourier space gives a direct pointwise estimate in physical space via ||(/'||c"' < C'||(/>||£i. 
On the other hand, working in Hilbert spaces H'^ is conceptually simpler since it allows for going 
back and forth between physical space (for the nonlinearity) and Bloch space (for the linear part). 
Moreover, localization of the resulting gap solitons follows directly in spaces. By the embedding 
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\(jk < C||0||//3 for k < s — 1 pointwise estimates are still obtainable although these are typically far 
from optimal. 

4.1 Preliminaries 

We have the asymptotic distribution 

Cin < ojn{k) < C2n VnGN,VA:GT2 (4.1) 

of bands ujn{k), with some constants Ci, C2 independent of k and n. This follows from the asymptotic 
"density of states" (see p. 55 of [22j) A^(A; k) = aX + 0{X^) as A — )• 00, where A^(A; k) is the number of 
eigenvalues Un{k) smaller than or equal to A. 
We introduce the diagonalization operator 

V{k)k(zj2 : 4){k]x)^4){k), ^„(A:) = {4){k--),pn{k;-)^^ 



Based on (4.1) we may estimate P. Similarly to |8i Lemma 3.3] we find that T>{k) for all k is an 
isomorphism between i7''(P^) and s > 0, where 

:= {J:||?|||:=5]|0,|2(l+jr <oo}. 

Moreover, ||P(A;)||, ||P(A;)~^|| < C independent of k. Thus, 1— t- </> is an isomorphism between 
L'^{T'^,H%F^)) and L'^iT^J^) and therefore we have 

Lemma 4.1 For s > the map 1— )• = VTcj) is an isomorphism between ff'^(M^) and := 
L^(T^,£^), = fj2 J2n<^N \'l^n{k)\'^ {1 + n)^ dk , i.e., we have the equivalence of norms 

ciUWhs < \mxs<c2\mHs. (4.2) 

For the sake of convenience, and since the idea of the proof is needed below, we also include a simple 
version of a Sobolev embedding theorem. 

Lemma 4.2 For </>,V G H%R'^), s > d/2, we have </',V' G , UWc^ < C'll'/'lk-, and < 
CUWhsUWhs. 

Proof. Let w{k) = (1 + \k\Y. Then € for s > d/2 and thus 

UWco < CUhl = C\\w-'^w\\l^ < CWw-'Wl^UWl^^ < CUh^^. (4.3) 
Next w{k) < w{k — i) + w{£) and thus 

\w{k){(i>*i>)ik)\ < [ \w{k-i)4>{k-i)i>{i)\di+ [ \w{£)(f>{k-i)^P{i)\di. (4.4) 
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Therefore, using Young's inequality ||(/) * < C'II'AIIlp IIV'IIl'?) 1/p + 1/'7 = 1/?' + 1, with r = 2,p 
2,q = 1 we have 



<CU*^l^\\^<C(^\\\w(^\*m\L2 + \M*\wi;\\\L2^ 

< c (^\\w4>h4i'\\L^ + Wwi^WML^) < 2C\\4>h^jmL2 < c\\nH^\m\H^- 



(4.5) 



□ 



By Lebesgue dominated convergence, <j) £ with s > 1 also implies (/'(x) — t- pointwise as 
|x| — ?■ oo, which shows that solutions in are indeed localized. Moreover, using the product rule, 



Lemma 4.2 can be generalized to ||i?^||(7fe < C||(/)||//s, k G N, k < s — d/2. Finally, by Lemma 4.1 these 
statements directly transfer to X'^. 

Lemma 4.3 Let s > 1 and e . Then (f>ip G X"^ , 4> G C^(M^) for k < s - 1, and 4>{x) — > as 
|x| — >• oo. 

To justify Taylor expansions of the spectral bands and the Bloch function we use the following 
lemma, which is proved in [39j. 

Lemma 4.4 Let V{x) G L^(P^), then for each n G N the band a;„(/c) and the Bloch function Un{k; x) 
are analytic in k on A; G \ Z„, where Zn is the set, where the algebraic multiplicity of u)n{k) as an 



eigenvalue of (2.1) is higher than one. 



4.2 Justification Step I: The extended Coupled Mode Equations 



To justify the general stationary CMEs (3.32) as an asymptotic model for stationary gap solitons near 



an edge a; = w*, we again consider (3.17), i.e. 

L{k;x) — uj^: — e^Q. (j){k;x) = —a (^(p * g (j) * g (pj {k; x) . (4.6) 

In contrast to the formal derivation of the CME in ^we now want to keep track of higher order remain- 
ders. We first apply the diagonalization operator V : (p{k;x) — )• (j){k), (pnik) = {4>{k;-),pn{k,-)) i2(j2~^ 
to get 



[u]n{k) - - e^^] (pnik) = -agn{k) 



(4.7) 



with gn{k) = {{(f> *B (j) *B 4)){k; ■),Pn{k] ■)) 



Lemma |43| allows us to consider (4.7) in the space X^ , s > 1. The multiplication operator 0Jn{k)—uj^ 



is not invertible since it vanishes at N points {n,k) = {ni,k^^^), . . . ^injq ,k^^^). This dictates our 
generalized Lyapunov-Schmidt decomposition of the solution (p, namely 



N 



^{k) = e-i^LS (A:) + V(fc), where e-^ff^^{k) = e'^ ^ ^ 



k + m — k^^^ 



(4.8) 



j=l meMj 
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an ansatz of the form (3.13) (with XDj{k + m)Aj ^ 

k = k^^\ since in e~^'fj^^{k; x) we have pnj{k; •) instead of pnj{k^^\ •) 
However, we have 



replaced by Bj [ k+^^k'^'^ 



with supp Bj C -D^r-i , ijjrij {k) = for /c G := U{Di : / G {1, . . . , N},ni = nj}, wihere Di is defined 
in (3.12) and where e„^. is the unit vector in the nj direction in M^. 

Remark 4.5 Note that in general £~^^^{k) is not the diagonalization of the leading order term in 

), except at 

(4.9) 
(4.10) 



with 



£ ^V'{!s(A;,x) = £ 



N / 

-^j;p„^.(fc(^),x) bA 



k + m — k^^^ 



and where for Bj G with s > 1 



N 

l|p||L2(ir2,_ffs(p2)) < C£ ||-Bj||j,2. 

T/izs follows from writing k = k^^^ — m + {k + m — k^^^), expanding 

PnAk,x) = PnAk'-^^ -m,x) + VkPnAki^\x) ■ {k + m- k^^'') 



(4.11) 



withk^^] G [min(/i;p^ — m^, A;/), max(/cp^ — m;, /c;)], I = 1,2, and using pn^ik^-'^ —m, x) = pnj{k^-^\x)e^^ 
which yields 



Si) 



N 



p{k,x)=£-'Yl E 4- (fc + m - fc(^')) • VkPn,{ki^\x). 

j=l m£Mj 



To prove (4-11), we may fix some (of the finitely many) j,m. Since uJnj{k) are simple for k G Dj, we 



have sup^g^. \\VkPnj{k, ■)\\h^{¥'^) ^ C, and it remains to estimate 



L2(T2) 



B 



k+m-k^J^ 



\k + m - k'~^\'^ dk 



<C£^ I \Bj{K)f\K\^ dK < C£^\\Bj\\l2 



-2|| 5 Il2 



E'GM2 



(4.12) 



for s > 1, and this proves (4-11)- 
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With the decomposition (4.8) we obtain the Lyapunov-Schmidt equations 

-(^XoXk)gnj{k), j = l,...,N, meMj, 



-{ujn^ik) -uj^- e u)Bj I 



N 



(4.13) 

(4.14) 



The goal is to solve (4.14) for the correction -0 as a function of B = {Bj)^^^ G L (Dgi—i, C ) and plug 



this into (4.13). It turns out that the right norm for Bj is ||2,2(£) where we recall 

ii + \v\YB 



B 



'j\\L2{D,r-l)- 



Note that L'^{D^r-i) = Ll{D^r-i) as spaces for any s > 0, but below we need the estimate < 
C||Sj||/^2(£) with C independent of e, cf. (1.5). 



Lemma 4.6 Let s > 1 and V G i?,^ ^^'^(M^), 5 > 0. For B G Ll{D^r-i) and ^ E s > 1, we have 



\~9\\xs<C 



loc 

N 



N 



AT" 



(4.15) 



Proof. Relying on the norm equivalence in (1.5), we derive most of the estimates in physical variables. 
Due to (4.8) and (4.9) we have g = \(p\'^4' = \^~'^'>PhS + /° + V'P('^~^'0ls + P + V')- We need to estimate 



norms of terms of the types 



(^^VS)'/, pV, pV-', ^nd^P^. 



(4.16) 



First note that e ^xp^^{x) = e ^^-^ i?j(ex)ii„^. (/c*^-'^; x). Below we implicitly use ||i?j(e-)un^ (fc*^-'); •)||hs < 



. {k^^^; •)|lcr''l ll-Si(£")ll-H''', which holds by interpolation, see, e.g., [351 §4-2]. Next, ||m„^ {k'^^'>; •)|| 



< 



j^fsi+i+i < C||F|| fsi-i+i for all 5 > 0, where the first inequality holds by Sobolev em- 

loc 

bedding and the second one by the differential equation. 



In estimating all except the first term in (4.16) we use the Banach algebra property in Lemma 4.2 
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For the first two terms we need to estimate ||(ei?j(e-))"'||//s for n = 2,3. We liave 



WieBjie-)) 



n||2 



< c 

< c 

< c 

< c 



{l + \k\r\T{ieB,{e-)r){k)\'dk 



\k\' 



dk 



\\{£Bj{e-))m2 + e'^'^-^+^' J \K\'^'\Bf{K)\'^dK 

^2(n-l)||^_^.||2(n-l) 
2(ra-l) 



2n-2+2s|| r)n||2 



L- \\^Bj{e-)\\^2 +e 

e2("-l)p.||5;r^)||i?,-||2,+e2(n-l)+2.p.„2n 



and hence 



||(eSj(e-))"lk= <Ce'^-'\\Bj\\'^s for n = 1,2, 3. 



(4.17) 



Note that for n > 2 this is much better than the naive estimate ||(e-Bj(e-))"||jys < C||ei?j(e-)||^s < 

(4.18) 



C||i?||^s based on (4.17) with n = 1. Next, for the third term in (4.16) we get 

\eB,{e-)f{-)\\Hs <Ce\\Bj\\Hs\\f\\Hs, 



and this together with (|4.17p can be used to prove (|4.15p. To show (|4.18p, we start with 

1 



\eBj{e.)fi-)\\Hs<\\eB,{e-)f{-)h^+C 



k 



-B, 



/(• 



L2 



The first term is estimated as \\eBj{e-)f{-)\\i^2 < ||oo ||/||l2 < ||/||j;^2, and for the second 

we note that w{k) < ew{^j^) +w{l) where w{k) = Thus, similarly to the proof of Lemma 
have, using Young's inequality. 



4.2 



we 



w{k){-^Bj[-)*f{.) 



< C 



<c 



w 



B, 



L2 



L2 



W 



B, 



Li + 



-B, 



1 



-B, 



*\Of{-)\ 



L2 



LI 



1^3 (i; 



< Ce\\B^\\L2, and H/H^i < CWfU^ (see (4.2)) yield 



Now (-) 13 j (-) 11^2 < e||4-|lL2, 
( 4l8l ). 

The 4th, 5th and 7th term in (4.16) are estimated simply using Lemma 4.2 and (4.11). The 6th 
term is treated the same way and is bounded by ^^(X^jLi ll-Sj||i;,2(^ r-i))^ t\i&i it is of higher order 
than the first term on the right hand side of (4.15). □ 



We seek now a small solution of (4.14). First, for e sufficiently small we have 



min \u}n{k) — > Ce 



2r 



(4.19) 



A;6supp(</)) 
neN 
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due to dkj^uJrijik^-'^) = dk^^Ujik^''^) = and the definiteness of the Hessian of Unj at k 



k^^\ We 



rewrite now (4.14) as 



i,n{k) = {ujn{k) - Lo^)-\-agn{k) + e'^ni^nik)) =: Fn^k) 



and solve ip = F{ip) on a neighborhood of 0, namely on B^v := {ip £ : < £^}j'n > 0, via the 

Banach fixed point theorem. 



Performing similar estimates as in the proof of Lemma (4.6) and using (4.19), we get 

|2 



H' 



+ 



where C = C(^^]^ ||i?j||j:/sm2'), \^\). The contraction property thus holds if 



+ \mH^)+e 



-2r 



for all '4^,4) & and if F maps Bg^v to Bi^v, i.e., using (4.15), if 



2 , .-2n,J,„3 ^ c^rj 



(4.20) 



(4.21) 



for all ip £ Bi;v. 



Condition (4.20) is satisfied if r < 1 and rj > max(2r — l,r) and (4.21) holds if max(2r — l,r) < 

(4.22) 



r] < 2 — 2r and r > 1. In combination these yield 



r/ G (r, 2 - 2r), < r < |. 



Here is the reason for the upper bound in (3.11) 



In conclusion we have the existence of ip satisfying 

AT 



< Ce^-^'^J^W^jWLUD^r-i) with < r < 2/3. 



(4.23) 



We now turn to (4.13). Plugging (4.8) into (4.7), truncating over k £ Dj and mapping k £ Dj — m 



to p £ D^r-i via p = e ^{k + m — k^^^) yields the so called extended CMEs (eCMEs) in the form 



QBi 



ldlu,n,{k^^^)pl + \dlu:n0'y)pl + 9fc,afc,a;„^.(fc(^))pip2 ) B, - Q, = e^R,{p), (4.24) 



j = 1, . . . ,N. Here Qj denotes the nonlinear term Mj in (3.32) with Ai replaced by Bi and truncated 
on p G -D^r-i , i.e., 



with ijj^^ given in (4.10). e^Rj{p) denotes the remainder, where we set r = min{3r — 1, 2 — 2r, 1}. The 



a 



,(p(^-'-))((V;S *B M *B As)(£p(^''") + k^'^ - m; ■),PnM'^-^ 
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remainder has the form 



(4.25) 



where q^^\p) is some cubic polynomial, D^^^LUn is a third derivative of w„. and ki''^ G is fixed. The 



first term in (4.25) comes from the Taylor expansion of ujn- and is bounded by 



3r-l| 



B 



enforcing r > 1/3 as set in (3.11 ). The leading order contribution to the remaining two terms in e^'Rj{p) 
comes from ijj, which is estimated via (4.23) to be ©(e^"^*"). This is guaranteed to be small since r < 2/3. 
The last contribution to e^Rj{p) comes from the ■^-independent terms in xd^,—i {p)9nj (ep + k^-'^ — m), 
i.e., from 



^XD^.-i (p) ((#s *s is *s is)(ep + k^^'^ 



m; ■),Pn, {ep + 



k^^^;-) 



e )L2( 



This difference is estimated in Lj to leading order via ||p||i:,2, which is 0{e) according to (4.11). 



Our first main result follows immediately from the above Lyapunov-Schmidt discussion. 

Theorem 4.7 Lets> 1,V G hI^}~^^\m?) for some 5 > and let ^ < r < There exist eo,Ci,C2> 
such that for all < e < Eq the following holds. Assume that there exists a solution (Bj)^^-^ G 



L^{D^r-i) of the extended CMEs (4.24) with \\Bj\\L2(^u^^_^) < Ci. Then ( [L2| ) has a solution (j) G H^{W) 
with 



N 



) - eY,B,{e-)un,{k^'^;-)\\Hsm <C2{e 



(4.26) 



where UnAk^~'^;x) are the resonant Bloch waves and Bj 



F-^B, 



Proof. By construction, a solution {Bj)jLi of (4.24) yields via (4 



estimate (4.23) and the norm equivalence in Lemma 



4.1 



because e^^^-i^Bj{ex)unj{k^-'^;x) =e ^^ls (^)' '^^^^ ^ "^%s (^) — ^ X^j=i (^'■"''*' ^) — Pi^)^ 



yield 



I a solution (j){k) of (4.6). The 

^-lJO)|| 



which is bounded in (4.11) and the approximation result (4.26) follows. 



< Ce^-^\ And 
x), 
□ 



4.3 Justification Step II: Persistence 



The remaining step is to make a connection between solitary waves of the CMEs (3.32) and the eCMEs 



(4.24). To obtain existence of solutions of the eCMEs (4.24), we show a persistence result of special 



so called reversible non-degenerate CME solitons to the eCME. This is quite similar to |131 §5] but in 
order to deal with an arbitrary set of extrema k^^\j G {1, A^}, the definition of reversibility needs 
to be generalized. 



The eCMEs (4.24) differ from the CMEs (3.32) in three ways: the Bj{p) are supported on D^r-i, 



the convolutions are truncated on D^r~i, and the remainder e^'Rj{p) is included. The idea is to handle 
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these differences as perturbations and thus seek a solution B = (-Bj)j=i,...,Ar of the eCMEs near a given 
solution A = (Aj)j=i,...,Ar of the CMEs. Note that the Aj are not supported on D^r-i and thus first 
must be truncated. 

We start with a formal discussion. We write the CME in abstract form as F(A) = and the eCME 

as 



XD,_,F(B) = /R(B), 



(4.27) 



assume a solution A G [Lq{W)]^ for any g' > of the CME, and look for solutions B G [Lj{D^r~i)] 
of the eCME in the form B = A^ + b with A^ = xd r-i-^ suppb C D^r-\. This yields 



N 



Jjb = N(b), p e with 

Je = XD._i^AF(A) and N(b) = e-R(A, + b) - (xd „_,F(A, + b) - J,b). 



(4.28) 



Note that in (4.28) we may replace N(b) by XD^r-i-'^(^) to display explicitly that p G D^r-i. 

We have F(Ae + b)-Jeb = F(Ae) + (L»^F(A£) - Je)b + G(b) with G(b) quadratic in b. Moreover, 



F(A,) = L^ME A - ._.AA(A) - {M{A,) - Xd ,._.AA(A)) = -(AA(A,) - xd ._,AA(A)), 



where Lp"^^ denotes the l inear part of the operator in (3.32) and denotes the A^— dimensional vector 



of the nonlinear terms in (3.32). XD r-i-^i-^)ip) a sum of convolutions Aj-^ * Aj^ * Aj^, and hence in 
AA(Ae) — Xd ,._iA/'(A) this yields terms of the form 



Pi ■J\P2\>£^ 



XD^r-1 (P - Pl)XD^r-i (Pl - P2)Ah (P - Pl)^h (Pl - P2)Aj3 iP2) dpi dp2 



which can be bounded by Ce"^ for any g > in Ll{D^r-i) due to the fast decay of A. Therefore, 



|F(A,)||^2 < Ce", and || (i?, F(A,) - J,)b||i2 < Ce'?||b||^2 



(4.29) 



by a similar estimate. 

Thus 1 1 N (b) 1 1 < C(e'' + e'^ + e'" 1 1 b 1 1 ^2 + e"? 1 1 b 1 1 + 1 1 b 1 1 ^2 ) ) and this suggests applying the contraction 



mapping theorem to (4.28) in the form 



b = J-iN(b). 



(4.30) 



To discuss J~^, we start with J : H'^{R'^) — )- H'^~'^{R'^). The continuous spectrum (Tc(J) of J equals that 
of L^'^^. Thus, if w = a;=K+e^f) is in a gap, then i7 and the quadratic forms defined by (A;(-^^)pf + 

^dl^uJnjik^-'^P2 + dkidk2(^nj{k^-'^)piP2, j=l> ■ ■ -jN have opposite signs such that a^J) is bounded away 
from zero. However, the problem is that J has a nontrivial kernel since Ker J contains at least 
dy^A-^dy^A. and iA which follows from the translational and phase invariances of the CME. For J^^ : 
Ll{D^r-i) — L^^2(-^e'— i) (if it exists) this implies that it cannot be bounded independently of e. 



The solution is to consider (4.30) in a subspace X^^^ C L1{D ^r-i) where ^ is bounded, and where 
b G ^rcv implies N(b) G X^^w This can be achieved by symmetries of the problem (1.2) if we assume 
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that J on H^(K ) has only dy^A, dy^A. and iA in its kernel. 

The original problem (1.2) is equi variant under the symmetries 

(t){xi,X2) = il(l){-Xi,X2) = (,24>{xi, -X2) 



(4.31) 



where 6,6 = ±1, i.e., (6,6) = (1,1) or (6,6) = (1,-1) or (6,6) = (-1,1) or (6,6) = (-1,-1), 
and similarly 

4>{xi,X2) = 60(a;2, xi) = 6<^(-a;2, -xi), (4.32) 

where again 6,6 = =tl- 

These symmetries have their counterparts in symmetries of the CMEs. As we show, the reversibility 
in the following definition firstly guarantees that b G Arcv implies N(b) G Aicv, secondly provides 
a leading order approximation erj^^^ satisfying (4.31) or (4.32), and lastly (under a non-degeneracy 
condition) guarantees the existence of a Gross-Pitaevskii solution cj) satisfying ( |4.31 ) or (4.32). 



Definition 4.8 A solution A of (3.33) is called reversible if it fulfills one of the symmetries in (4.33) 
or ( |4.34D . 

^i(y) = si^i"(— 2/1, 2/2) = S2^i'(?/i, —112) for all i G {1, . . . , A^} with si^2 = ±1 independent of i, 

(4.33) 

and where the indices i' and i" are defined by 



k^^') = {-kxKkf ) if kX' <\ and i' = i if k\ 
k^^") = {kf> ,-kf ) if kf> <\ and i" = i if k. 



1 

2' 



Ji) _ 1 
" 2- 



Ai{y) = Si^j'(y2,yi) = •S2^j"(— ?/2,— yi) for all i G {l,. . . , A^} with si^2 = il independent of i, 

(4.34) 

and where the indices i' and i" are defined by 
. k^) = {k^\k^), 



fc(^") = (_fc«,_A;f)) 



^/ ^1,2 < 2' 



i" = i 



if A:W = M, 



= (_A;«,A;P) ^/ A:« = i,A;«/M, k^'") = {kf , -kf) if kf = \,k^>^^M. 
We define the space 

Xrev = {f G lI{D ^r-i) : f satisfies (Osl) or (Ol}. 



Definition 4.9 A solution A of (3.33) is called non- degenerate if Ker J = {9^^^ A, A, i A}. 
We explain next that under the assumption A G Xrev we, indeed, have the property 



b G Xr, 



N(b) G Xr. 



(4.35) 



Note that this resembles the question of inheritance of symmetries of the full problem for (f) by the 
Lyapunov-Schmidt reduction, see e.g. Prop. 3.3 in [18]. We are, however, using a generalized Lyapunov- 
Schmidt reduction in which a whole neighborhood of the kernel is isolated. Secondly, we wish to carry 
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symmetries of the scalar problem for <j) over to symmetries of the vector problem for the envelopes B. 



We inspect, therefore, problem (4.35) by hand. For the sake of brevity we present the analysis only for 



the first symmetry in (4.33), i.e., 

Ai{y) = siAi„{-yi,y2) for alH G {!,..., iV}. 



(4.36) 



This is equivalent to Ai{p) = siAi"{—pi,p2) for all i G {1, . . . , A^} in Fourier variables. If we denote 



SA{y) = ib(yli//, A2", . . . , ANri){—yi, 7/2), condition (4.36) reads A = SA. 

In the term XD r-i'^i-^e + b) in N(b) the operator F is the CME operator, which commutes with 



S. In the argument of F the function satisfies (4.36) because A does and the symmetric cut-off 
-A-e = xd r-i (p)-^ preserves the symmetry. The linear operator commutes with S because it is 
a symmetric cut-off of the linear operator corresponding to CMEs. It remains to discuss the term 



e'"R(Ae + b) in N(b). The first term on the right hand side of (4.25) clearly commutes with S and 
so does the last one because it is a symmetric cut-off of the nonlinear part of the CME operator. 
The term gnj{k) is more complicated because it involves pnj{k;x) in the convolutions instead of the 
A;— independent pnj{k^-'^;x). After solving (4.14) for the regular part -0 in terms of B, we obtain that 
Qn^ (k) consists of the terms 



fais^j^ XD^{k + mmo.*BCl3*Be^)ik;-),Pn,{k;-))L2(^p2^, 



(4.37) 



where = Pn^(.k;x) Ylm&M^ XdA^ + m)Ba 



and = Pn^{k;x)Y,meM^X^DAk - m) 



k-m+k'-°''l 



We have 



m£Mj n,o,q£Aa,i3,-,,j,m 



^ n nam V / 



with 



w. 



noqm 



(k) =XDAk + m) / / XD^{k-r+n)Ba 



k-r+n-k("'> 



XDgir-s+o)Bij 



r—s+o—k^^^ 



T2 T2 



X X-D^{s-q)B^ ( ) {pnAk - r; ■)Pnp{r - s; ■)pn.,{s; ■),Pn,{k; •))l2(p2) dsdr. 



Employing the same change of variables as in (3.24) and using (3.29), we get 



fj-naqmip) ■= wfoqli^P + ^^^^ - m) = XD^r-i (p) J J XD^r-i iP " f)Ba{p - r)XD^r-i {r - s)Bp{f - s) X 

X XD^.-i (5)^7(5) (pn,(/c(°) -n + e{p-f);-)pnf,{k^^^^ - o + e(f - s); •) x 
X ^(-fc^^-* + q + es; ■),Pnj{ep + k^^'> - m; •)>L2(p2) dsdf, 

(4.38) 



where = D^.r-. n T^-fc<^'+fcW+o-, ^ ^ i^+k^ 



For the symmetry S we need to consider fJ.noqm{—pi,P2) = fj-noqmipi, —P2)- Assuming that b 
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(and A) satisfies ( |4.36 ), we get that also B does. Thus, using Bi ((-pj) ~ ^) ~ {{ ~^ ^) 
siI3iii (p + ( 1}^ )) and the change of variables f — >• (fi, — f2), s — )■ (si, — S2), we get 



x(lv:(^(")-n + e((-J-(4));-)p;^(fc('^)-o + .((4)-(4));^^ 



+ '7 + e(4);- 



L2(P2) 



dsdf, 



(4.39) 



where = L>2e'-i n e'^ 



J/3) 



(/3) 



1X1) 



+ ( ) + ( 4 ) - ( 4 ) ) ' ^ = n e-i (T2+ 



91 ■ 
-92 , 



Finally, we use the Bloch function symmetry Pn((A;i, — ^2); 2;) = pn{k;{2TT 



xi,X2))- In our case this means 



(a) 



(_"4)+.(p-f);(2--0 • 



■y" = n„ 



If ^2 £ (~ 1/2, 1/2), then n2 = and we have for k^°' ^ = ( ^^^^ J and n^'^ 

( -4 ) + ^(P - r); ( ^-/i )) = Pn^„ - n + e(p - f ); ( ^--^ )) . 

If — 1/2, then 722 G {0, 1} and ( ^^^j J — ( .I^^^ ) = k^°"> — h for some n G with n 7^ n. Thus 

- ( _^ ) + eip- r); ( ^--^ )) = l^n. (a^^^^ - n + e(p - f ); ( ^--^ )) 



Pn. 



so that a" = a. 



Performing the same analysis for the Bloch functions Pnp^Pn^ and pn^ in (4.39) and using sf = si, 
we get 



Anogm(-Pl,P2) = ^1 A^^^"^" (p) , 

where some of the doubly primed and 'tilded' indices may equal the bare ones and where n G Mq,, o G 
M^, q G M^, and m G Mj. After the sum in n, o, g and m, we get 



r^^^'(-l^l,P2) = Sl/""''"^"^" 



In conclusion N = SN and (4.35) holds. 

With Ker J = {dy^A^dy^K^iA^ the operator has an 0{1) bounded inverse on Arev, i-e., with a 
bound independent of e, and ( 4.30| ) can now be solved producing 
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Theorem 4.10 Let s > 1. There exists an eq > such that for all < e < Eq the following holds. Let 
Lo = io^: + e'^Q be in a band gap, let A be a reversible non- degenerate solution of the CME (3.33) with 
A G L^(R^) for all q > 0, and let 1/3 < r < 2/3. Then there exists a C > and a solution B of the 
eCME such that 



IB- A| 



i) < f = min{3r - 1, 2 - 2r}. 



(4.40) 



Corollary 4.11 The solution (j) constructed in Theorems J^.l and J^.IO is a localized solution of (1.2), 
it is symmetric according to (4-31) or (4-32), and 



N 



EY,M^->n,{k^'^r)\\H^ 



< Ce'. 



(4.41) 



Proof. We first show that the reversibility of B and the symmetry of Bloch functions provides a e ^??ls' 
that satisfies (4.31) or (4.32). Let us work out explicitly only the first symmetry in (4.31). Recall that 



^LS= / e^'-^X;Pn,(A:;x) XD,{k + m)B 

/T2 

^ j=l mSMj 



dk. 



Using pn^{k;{-xi,X2)) = Pn^„{{ki, -k2); x) with nj» = nj, Bj{p) = siBj//{-pi,p2) = siBj//{pi, -p2) 
together with the definition k^^"^ = {k[^\ —k2'') if k2^ < ^ and k^^"^ = k^^^ if k2^ = ^, we get for 



uU) ^ 1 



N 



V^!i{-xi,X2) = si / eHfci.-fe)- Vlv7((^i,-^2);x) V XD,{k + m)B, 

/T2 — — 

^ j=l mGMj 



dk. 



Because m2 = Vm £ Mj when A;2"'^ < ^ and since Mj = Mjn, and XDj{k + rn) = xd 
we have after the change of variables k — t- (ki, —k2) 



ki 

j" V V -k2 



m 



N 



Vis{-xi,X2) = si I ^e'''-^Ypn^„{k;x) ^ XDy,{k + m)Bj 

j=i 



J2 



k+m—k^^ ' 



dk = r/[°s'(x) 



For ^2 = ^ we have 



N 



rffl{-xi,X2) = si / e^C^i' '=2)-^Vpn ((/ci,-A;2);a;) V XDj{k+m)Bf' 
Jt2 



dfe 



because ^ '^^(^^ ^ = A;^-' ^ — ( 5 ). Next, using firstly ( _^^2 ) + ( i ) = with some m = Mj,m / m, and 



33 



secondly xz), {k + m) 



XD^„ 



-k2 



N 



+ m ) , we arrive (after the change of variables k — )• {ki, —k2)) at 



J2 



Next, V in (4.14) inherits the symmetry, so that cf) is symmetric. The estimate (4.41) follows from 

□ 



the triangle inequality with (4.26) and (4.40) 



Remark 4.12 In (4.40) we use f = min{3r — 1, 2 — 2r} although below ( 4.24[ ) we defined f := min{3r — 
1, 2 — 2r, 1}. This is because for 1/3 < r < 2/3 we have min{3r — 1, 2 — 2r, 1} = min{3r — 1, 2 — 2r}. 

The optimal value of f is f = 4/5 attained at r = 3/5. Based on the formal asymptotic expansion 
in (3.13) and (3.20), we see that the next order term -i^^-*^) (just like 'ijj^'^^) consists of terms of the type 
{k^^^;x), where F is an envelope and q a periodic carrier wave. therefore, consists of 

= 0{e). As a result, the formal asymptotics predict 



F 



■,x e 



and 



(1)1 



terms e'^F{ex)q{k^^^ 

convergence. Thus, while the estimate (4.41) guarantees convergence of the CME approximation. 



it does not appear to be sharp. But if all third derivatives of conj vanish at k^^\ like for separable 
potentials [13], we have f = min(4r — 1, 2 — 2r, 1) with 1/4 < r < 2/3 and the optimal value is f = 1 
attained at r = 1/2. It is, however, unclear which non-separable potentials result in vanishing third 
derivatives of the bands at gap edge extrema. J 



Remark 4.13 As said in the previous remark, the formal asymptotics predict that the error ip in ( |4.8| ) 
has the form 

ip{x) = e^F{ex)w{x) 

with F e i?'?(R2) for all g > and w{x) e cl^'^i^^) with w{2'k,X2) = e^'^'^^w{Q,X2), w{xi,2'k) = 
e^'^^^'^w{xi,Q). In this case 

Hs = O(e^) implies ||V||l- = 0{e^+^). 



loo > Cie° with Q < 1 + i.e., iV'la^o)! > C'le" for some xq G M^. Since ip 



x)\ > in some neigborhood of xq of diameter 6i > 0, say. Moreover, 



To see this, assume that 
is continuous, this implies 
since F is continuous, |e^F(ex)| > for |x| < 62/£- Therefore, for e sufficiently small, ^ 
for X in 0{5i) wide neighborhoods of Xj = xq + (2ji7r, 2j27r), where j'l, j2 G Z and |ji|, \j2\ < C2/e, see 
Figjrjfor a ID sketch. Since in 2D there are at least Cse"^ such neighborhoods, we obtain 



(4.42) 



which contradicts HV'llii'^ ^ Ce^ as e — )• due to a — 1 < ^. 

To make this argument rigorous for the error ip = (j) — e^^^^i Aj{e-)unj{k^-'^; •), we could split off 
the next term in the formal asymptotic expansion and show that the remainder is of higher order, i.e., 
■0 = e^V^^"* Then we can estimate = 0{e^~'^^') using the analysis from ^4.2, but here 

we refrain from these tedious calculations. J 
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1/e 




1/8 many 0(1) intervals where |\|/|> C e" 



Figure 7: Sketch of the argument for HV'IIh- = 0{ef^) = ©(e^+i). 



5 Numerical Results on Reversible Gap Solitons 

We numericahy compute some representative cases of gap sohtons and their asymptotic envelope ap- 
proximations e(/>*^°^(x) = eY^^^^Aj{ex)unj{k^-'^]x). We do not attempt to provide an exhaustive study 
of possible GS solutions but rather select only several cases to corroborate our analysis. Namely, we 
select GSs bifurcating from the edges S2 and S5. The latter case is of particular interest as it fea- 
tures a situation whose occurrence is impossible for separable potentials V{x). To our knowledge this 
case has not been studied before and the presented GSs are novel. We also check the reversibility 



and non-degeneracy conditions which are sufficient for persistence, see ^4.3, In addition, we com 



pute the convergence rate in e, i.e., in the square root of the distance to the gap edge, of the error 



6g^s--e</>(o)||H2. 



A 4th order centered finite difference discretization is used for (1.2). The computational domain is 
a square x E [-LIgs/2,I?gs/2]^ selected large enough so that the asymptotic approximation ecj)^^^ (x) of 
the GS is well-decayed at the boundary and zero Dirichlet boundary conditions are then used. Equation 



(1.2) is then solved via Newton's iteration using ecp^ ' as the initial guess. The computational domain 



is in practice reduced to its quarter using the corresponding reversibility symmetry. 



5.1 Gap Solitons near lo = S2 

Near the edge uj=S2 we limit our attention to real, even GSs and to symmetric vortices of charge 1. As 
the coupled mode system near ci; = S2 is a scalar nonlinear Schrodinger equation, see §3.2.2 



one can 



search for solutions of the form A{y) = i?(r)e'™'^, where r = -^\/y\ + 1/2 ' ^ ~ arg(yi+iy2), and m G N. 
We choose i? > and m = corresponding to the so called Townes soliton, and m = 1 corresponding 
to a vortex of charge 1. The function i?(r) satisfies the ODE 

i?" + -i?' + Oi? - ^i? - (T7i?3 = 0, (5.1) 
where i?(0) > 0, i?'(0) = for m = and i?(0) = 0, i?'(0) > for m = 1. For m 7^ the initial-value 



problem for the ODE (5.1) is ill-posed but can be turned into a well-posed one via the transformation 

Q" + ?^^Q' + l^Q - a7r2-Q3 ^ q (5.2) 



Q = r "^R{r) leading to 
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~ r'™! as r 



with Q(0) > 0, such that R{r 
method searching for Q{r) vanishing as r — )• oo. 

For m = we have the reversibihty A{—yi,y2 
with si = S2 = 1 since A is real. The non-degeneracy condition on J in Theorem |4.10| is known to be 



0. We solve equation (5.2) numerically via a shooting 
A{yi,—y2) = A{y), which is the same as (4.33) 



satisfied by the positive ground state A [2H |9] and conditions of this theorem are, therefore, satisfied. 

Figure [s] shows the profiles of the envelope A, of the asymptotic approximation (x) and of 
the GS (f>{x) computed via Newton's iteration on ( |1.2[ ). A GS deep inside the gap {s2,ss) obtained 
via a homotopy continuation in w from the (p{x) in Fig. [8] is plotted in Fig. |9]^a), while (b) shows the 
e-convergence of the approximation error. Here the e^'^^ convergence rate is better than the estimate 
proved in Corollary 4.11 and even better than the rate predicted by formal asymptotics in Rem. 

Km 
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Figure 8: Profiles of the even real GS at a; = S2 + e'^Q,, e = 0.1,0, = 1 . (a) A{y); (b) the corresponding 
leading-order GS approximation eA{y)ui{M; x); (c) the numerically computed GS at w = S2 + e'^Q. 
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Figure 9: (a) Profile of a GS corresponding to the even real family that bifurcates from a; = S2 in 
Fig. [sj The plotted GS is deep inside the gap (s2,'S3) at a; 1.78 (corresponding to e 0.28). (b) 
e-convergence of the error \\(\) — e^'^^^\ii2( 



For m = 1 the solution is complex and we have the reversibility A(—yi, 7/2) = —A{yi, — ^2) = ~A{y), 
which is (4.33) with si = — §2 = —1- Figure [To] shows the modulus and phase of the envelope A, of the 
asymptotic approximation e(f)^'^\x) = eA{ex)ui{M; x) and of the computed GS. The non-degeneracy 
of the envelope is illustrated in Fig. fllla), which plots the 4 smallest eigenvalues (in modulus) of the 
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Jacobian operator J of the CMEs evaluated at the vortex A: 3 eigenvalues converge to zero as the 
computational domain size grows while the fourth one stays bounded away from zero. Figure [ITj^b) 
presents the e-convergence of the approximation error ||</> — e(f)^^^\jj2(^-^2y The resulting convergence is 
very close to e^, which is the prediction based on formal asymptotics. 




Figure 10: Profiles of the vortex GS at a; = §2 + e^^^, e = 0.09, Q = I . (a) and (b) modulus and phase 
of A(y) resp.; (c) and (d) modulus and phase of the corresponding leading-order GS approximation 
eA{y)ui{M; x) resp.; (e) and (f) modulus and phase of the numerically computed GS at to = S2 + e'^^ 
resp. 
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Figure 11: (a) The four smallest eigenvalues of the Jacobian J in Theorem 4.10 at the solution A 
(a-d) for a range of sizes of the computational domain. 
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(b) e-convergence of the error 
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5.2 Gap Solitons near lo = 

We limit our attention here to gap solitons with real positive envelopes satisfying the symmetries 
Ai = A3,A2 = A4 and Ai{-yi,y2) = ^1(^1,-^2) = Ai{-y2,yi) = ^2(2/1,2/2), which is (|4.34|) with 



•si = S2 = 1 for each Aj. Such solutions of the CME system (3.39) can be found by first setting 



02 = and computing radially symmetric positive solutions Ai = A2 = A3 = A4^ = R{r), where 
r = \/ y\ + 2/2 ' ^ shooting method and then performing a homotopy continuation in 02 on 
the system of the first two equations in (3.39) employing the symmetry A\ = A3, A2 = A4 up to the 
original value 02 = 0.096394. 

We normalize the Bloch functions vi{x) := UQ{{kc,kc);x),V2{x) := UQ{{—kc,kc);x), V3{x) := 
U6{i—kc,—kc)',x) and U4(x) := UQ{{kc, —kc); x) so that 



t>2(-rEi, 2:2) = t'i(2;i,X2), V3{xi,-X2) = V2{xi,X2) and V4^{-xi,X2) = V3{xi,X2) (sec P). (5.3) 



This normalization implies that e^^*^) (x) is real and even in both variables. These symmetries are used 
to reduce the computational domain to one quadrant and restrict to the real arithmetic. 

shows the envelope Ai{y), the GS approximation e(j)^^^ and the computed GS (j). The 
is not radially symmetric due to the mixed derivative dy^dy^ in (3.39), 



Figure 



12 



12 



envelope Ai{y) in Fig. 

but looks radially symmetric because the coefficient 02 is relatively small (02 ~ 0.0964). Profiles of 
A2, ■ ■ ■ , A4 are not plotted as they can be obtained from Ai via the above mentioned symmetries. 
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Figure 12: Profiles of the even, real GS at = 55 + e^ri,e = 0.1,17 = —1 . (a) Ai{y); (b) the 
corresponding leading-order GS approximation sY^'^^^ Aj{y)vj{x); (c) the numerically computed GS 
at w = S5 + e'^Q. 

A closer look at the structure of (j) near the origin, an illustration of the non-degeneracy of A, and 
the e-convergence of the approximation error are provided in Fig. 13 The obtained rate is about e"'^^. 



which is once again close to the rate predicted by the formal asymptotics. 

6 Conclusions 

We have derived systems of Coupled Mode Equations (CME) which approximate stationary gap solitons 
(GSs) of the 2D periodic Nonlinear Schrodinger Equation/ Gross Pitaevskii equation near a band edge. 
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Figure 13: (a) Detail of the profile in Fig. [12] (c). (b) the non-degeneracy of A. (c) e-convergence of 
the error ||(/) - e(/)(°)||^2(R2). 

In contrast to [13] we make no assumption on the form of the periodic potential V{x). While in the 
case of a separable V{x) |13j the derivation is possible in physical variables, here in general it has to 
be performed in Bloch variables. We have rigorously proved via the Lyapunov-Schmidt reduction that 
reversible non-degenerate solitons of the CME yield GSs of the Gross-Pitaevskii equation. We have 
also provided an H'^{M?), s > 1 estimate on the approximation error showing that it is 0{e^^^) for 
GSs with the spectral parameter 0{e^) close to the band edge. The analysis has been corroborated by 
numerical examples including one which features novel GSs bifurcating from a band edge Bloch wave 
located outside the set of vertices of the first Brillouin zone, which is impossible in the case of separable 
potentials. 
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